//
// Grv98Pdf.cxx
//
// Created by TB on 11/23/11.
// Copyright 2011 BNL. All rights reserved.
//

#include <cstdlib> // For abs(int)
#include <sstream>
#include <stdexcept>
#include <vector>

#include "tppmc/Enumerations.h" // For LHAPDF::Flavour
#include "tppmc/Grv98Pdf.h"
#include "tppmc/Particle.h"
#include "tppmc/PartonContent.h"

// Access to FORTRAN GRV98PA subroutine (see grv98.f)
extern "C" {
   void grv98pa_(int* ISET, double* x, double* Q2, double* UV, double* DV,
                 double* US, double* DS, double* SS, double* GL);
}

namespace tppmc {

   Grv98Pdf::Grv98Pdf(int set)
   : mSet(set) {
      // Throw an exception if the set is not valid
      if(set < 1 or set > 3) {
         std::ostringstream oss;
         oss << "Grv98Pdf::Grv98Pdf(): " << set << " is not a valid set";
         throw std::invalid_argument(oss.str());
      } // if
   }

   void Grv98Pdf::Evaluate(const pythpart& parton, double scale,
                           PartonContent& partcontent) const {
      PARTONCONTENT p;
      Evaluate(parton, scale, p);
      
      // Intialise to zeros:
      partcontent = PartonContent(0.);
      
      partcontent[LHAPDF::UP]      = p.up; 
      partcontent[LHAPDF::DOWN]    = p.down;
      partcontent[LHAPDF::UBAR]    = p.anti_up;
      partcontent[LHAPDF::DBAR]    = p.anti_down;
      partcontent[LHAPDF::STRANGE] = p.strange;
      partcontent[LHAPDF::SBAR]    = p.anti_strange;
      partcontent[LHAPDF::GLUON]   = p.glu;
   }

   void Grv98Pdf::Evaluate(const pythpart& parton, double Q2,
                           PARTONCONTENT& partcontent) const {
      std::vector<double> dens(13, 0.);
      // grv98pa needs a non-const pointer to set:
      int ISET = mSet;//1; //LO
      double UV, DV, US, DS, SS, GL;
      double x = parton.X();
      grv98pa_(&ISET, &x, &Q2, &UV, &DV, &US, &DS, &SS, &GL);
      
      dens[8] = UV+US;
      dens[7] = DV+DS;
      dens[4] = US;
      dens[5] = DS;
      dens[9] = SS;
      dens[3] = SS;
      dens[10]= 0.;
      dens[2] = 0.;
      dens[11]= 0.;
      dens[1] = 0.;
      dens[12]= 0.;
      dens[0] = 0.;
      dens[6] = GL;
      
      //Assign parton densities PROTON:
      partcontent.up           =  dens[8] / x; //
      partcontent.down         =  dens[7] / x; //
      partcontent.anti_up      =  dens[4] / x; //
      partcontent.anti_down    =  dens[5] / x; //
      partcontent.strange      =  dens[9] / x; //
      partcontent.anti_strange =  dens[3] / x; //
      partcontent.charm        =  dens[10]/ x; //
      partcontent.anti_charm   =  dens[2] / x; //
      partcontent.bottom       =  dens[11]/ x; //
      partcontent.anti_bottom  =  dens[1] / x; //
      partcontent.top          =  dens[12]/ x; //
      partcontent.anti_top     =  dens[0] / x; //
      partcontent.glu          =  dens[6] / x; //
   }

   /**
    Returns true if the distribution of the input parton type
    is supported by this function.
    */
   bool Grv98Pdf::Supports(const pythpart& parton) const {
      return ::abs(parton.Pdg() < 4) or 21 == parton.Pdg();
   }

   std::string Grv98Pdf::Name() const {
      std::ostringstream oss;
      oss << "GRV 98 (" << mSet << ")";
      return oss.str();
   }
   
} // namespace tppmc
